library(tidyverse)
library(broom)
library(haven)
library(stargazer)

data <- read_dta("Study_1.dta")

# Create treatment variable
data <- data %>%
  mutate(treatment = case_when(
    trustExperimentSplit == 1 ~ 0,
    trustExperimentSplit == 2 ~ 1
  )) %>%
  mutate(treatment = factor(treatment, labels = c("Control", "Treatment")))

# Create trust variables
data <- data %>%
  mutate(
    trustMP = coalesce(trustExperimentGroup1MPRe, trustExperimentGroup2MPRe),
    trustMSP = coalesce(trustExperimentGroup1MSPRe, trustExperimentGroup2MSPRe),
    trustScotMin = coalesce(trustExperimentGroup1ScotMinRe, trustExperimentGroup2ScotMinRe),
    trustUKMin = coalesce(trustExperimentGroup1UKMinRe, trustExperimentGroup2UKMinRe),
    trustCivWest = coalesce(trustExperimentGroup1CivWestRe, trustExperimentGroup2CivWestRe),
    trustCivScot = coalesce(trustExperimentGroup1CivScotRe, trustExperimentGroup2CivScotRe)
  )

# Prepare data for plotting
plot_data <- data %>%
  pivot_longer(cols = c(trustMP, trustMSP, trustUKMin, trustScotMin, trustCivWest, trustCivScot),
               names_to = "outcome", values_to = "trust") %>%
  mutate(outcome = factor(outcome, 
                          levels = c("trustCivScot", "trustCivWest", "trustScotMin", "trustUKMin", "trustMSP", "trustMP"),
                          labels = c("Scottish Civil Servants", "UK Civil Servants", "Scottish Ministers", "UK Ministers", "MSPs", "MPs")))

# Calculate means and confidence intervals
means_data <- plot_data %>%
  group_by(outcome, treatment) %>%
  summarise(
    mean_trust = mean(trust, na.rm = TRUE),
    se_trust = sd(trust, na.rm = TRUE) / sqrt(n()),
    ci_lower = mean_trust - 1.96 * se_trust,
    ci_upper = mean_trust + 1.96 * se_trust,
    .groups = "drop"
  )

# Function to perform t-test and return p-value
t_test_p_value <- function(outcome_var, data) {
  t_test_result <- t.test(outcome_var ~ treatment, data = data)
  return(t_test_result$p.value)
}

# Calculate ATE and p-values
ate_data <- plot_data %>%
  group_by(outcome) %>%
  summarise(
    p_value = t_test_p_value(trust, cur_data()),
    .groups = "drop"
  ) %>%
  left_join(
    means_data %>%
      pivot_wider(id_cols = outcome, names_from = treatment, values_from = mean_trust),
    by = "outcome"
  ) %>%
  mutate(
    ate = Treatment - Control,
    stars = case_when(
      p_value < 0.001 ~ "***",
      p_value < 0.01 ~ "**",
      p_value < 0.05 ~ "*",
      TRUE ~ ""
    ),
    ate_label = sprintf("ATE = %.2f%s", ate, stars)
  )

# Create Figure 2
plot1 <- ggplot() +
  geom_point(data = plot_data, aes(x = outcome, y = trust, color = treatment), 
             position = position_jitterdodge(jitter.width = 0.3, jitter.height = 0.35, dodge.width = 0.8),
             alpha = 0.075, size = 1.5, shape = 16) +
  geom_pointrange(data = means_data, 
                  aes(x = outcome, y = mean_trust, ymin = ci_lower, ymax = ci_upper, group = treatment, color = treatment),
                  position = position_dodge(width = 0.8),
                  size = 0.8) +
  geom_text(data = ate_data, aes(x = outcome, y = 7.2, label = ate_label),
            vjust = -0.5, size = 3) +
  scale_color_manual(values = c("#11CE9E", "#CE1141FF"), name = "Group") +
  labs(y = "Trust score",
       caption = "* p<0.05, ** p<0.01, *** p<0.001") +
  theme_minimal() +
  theme(
    axis.text = element_text(size = 14),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    axis.text.x = element_text(angle = 0, hjust = 1),
    legend.position = "bottom"
  ) +
  coord_flip(ylim = c(0, 7.5))

ggsave("Figure_2a.png", plot1, width = 12, height = 6)

# Figure 6 (supplementary materials)
control_data <- plot_data %>% 
  filter(treatment == "Control") %>%
  mutate(outcome = factor(outcome, levels = c("MPs", "MSPs", "UK Ministers", "Scottish Ministers", "UK Civil Servants", "Scottish Civil Servants")))

control_dist_plot <- ggplot(control_data, aes(x = trust, y = ..prop.., group = outcome)) +
  geom_bar(fill = "#CE1141", alpha = 0.5, color = "#CE1141", size = 0.5) +
  facet_wrap(~outcome, ncol = 2) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(x = "Trust Score (Higher = More Trust)", y = "Percentage") +
  theme_minimal() +
  theme(legend.position = "none")

ggsave("Figure_6.png", 
       control_dist_plot, width = 10, height = 8)

library(kableExtra)
# Function to perform t-test and return tidy results with group means
run_t_test <- function(data, outcome_var) {
  t_test_result <- t.test(data[[outcome_var]] ~ data$treatment)
  control_mean <- mean(data[[outcome_var]][data$treatment == "Control"], na.rm = TRUE)
  treatment_mean <- mean(data[[outcome_var]][data$treatment == "Treatment"], na.rm = TRUE)
  n <- sum(!is.na(data[[outcome_var]]))
  
  tibble(
    outcome = outcome_var,
    mean_diff = t_test_result$estimate[2] - t_test_result$estimate[1],
    control_mean = control_mean,
    treatment_mean = treatment_mean,
    p_value = t_test_result$p.value,
    n = n
  )
}

# List of outcome variables
outcome_vars <- c("trustMP", "trustMSP", "trustUKMin", "trustScotMin", "trustCivWest", "trustCivScot")

# Run t-tests for all outcome variables
results <- map_df(outcome_vars, ~run_t_test(data, .x))

# Create a formatted table
table <- results %>%
  mutate(
    outcome = case_when(
      outcome == "trustMP" ~ "MPs",
      outcome == "trustMSP" ~ "MSPs",
      outcome == "trustUKMin" ~ "UK Ministers",
      outcome == "trustScotMin" ~ "Scottish Ministers",
      outcome == "trustCivWest" ~ "UK Civil Servants",
      outcome == "trustCivScot" ~ "Scottish Civil Servants"
    ),
    significance = case_when(
      p_value < 0.001 ~ "***",
      p_value < 0.01 ~ "**",
      p_value < 0.05 ~ "*",
      p_value < 0.1 ~ ".",
      TRUE ~ ""
    )
  ) %>%
  select(outcome, mean_diff, control_mean, treatment_mean, p_value, significance, n) %>%
  mutate(across(where(is.numeric), ~round(., 3)))

# Create and save LaTeX table
latex_table <- kable(table, format = "latex", booktabs = TRUE,
                     col.names = c("Outcome", "Mean Difference", "Control Mean", "Treatment Mean", "p-value", "", "N"),
                     caption = "Treatment Effects on Trust in Political Actors") %>%
  kable_styling(latex_options = c("striped", "hold_position"))

# Save the table to a file
writeLines(latex_table, "Table_2.tex")

# Print the table in the console
print(latex_table)